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Two-frequency radiative transfer (2f-RT) theory is developed for classical 
waves in random media. Depending on the ratio of the wavelength to the 
scale of medium fluctuation 2f-RT equation is either a Boltzmann-like integral 
equation with a complex-valued kernel or a Fokker- Planck- like differential 
equation with complex-valued coefficients in the phase space. The 2f-RT 
equation is used to estimate three physical parameters: the spatial spread, 
the coherence length and the coherence bandwidth (Thouless frequency). A 
closed form solution is given for the boundary layer behavior of geometrical 
radiative transfer and shows highly nontrivial dependence of mutual coher- 
ence on the spatial displacement and frequency difference. It is shown that 
the paraxial form of 2f-RT arises naturally in anisotropic media which fluc- 
tuate slowly in the longitudinal direction. © 2008 Optical Society of America 

OCIS codes: 030.5620, 290.4210 



1. Introduction 

Let Uj,j = 1, 2, be the random, scalar wave field of wavenumber kj,j = 1, 2, The mutual co- 
herence function and its cross-spectral version, known as the two-frequency mutual coherence 
function, defined by 

r 12(x , y) = ^ l( * + |_ ra * (1) 

where (■) stands for the ensemble averaging, is the central quantity of optical coherence 
theory, from which the two-space, two-time correlation function can be obtained via Fourier 
transform in frequency, and therefore plays a fundamental role in analyzing propagation of 
random pulses. 3 ' 4 ' 17,18 ' 21 The motivation for the scaling factors in ([1]) will be given below, 
cf. (H. 
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In this paper, we set out to analyze the two-frequency mutual coherence as function of 
the spatial displacement and frequency difference for classical waves in multiply scattering 
media. This problem has been extensively studied in the physics literature (see 2,17 ' 23,26 and 
references therein). Here we derive from the multscale expansion (MSE) the two-frequency 
version of the radiative transfer equation which is then used to estimate qualitatively the 
three physical parameters: the spatial and spatial frequency spreads, and the coherence 
bandwidth, also known as the Thouless frequency in condensed matter physics. Moreover, 
we show that the boundary layer behavior of the two-frequency radiative transfer (2f-RT) 
equation is analytically solvable in geometrical optics. The closed form solution (f43|) provides 
detailed information of the two-frequency mutual coherence beyond the current physical 
picture 23 ' 24 ' 26 (see the discussion about 

To this end, we introduce the two-frequency Wigner distribution whose ensemble aver- 
age is equivalent to the two-frequency mutual coherence and is a natural extension of the 
standard Wigner distribution widely used in optics. 6 ' 14 A different version of two-frequency 
Wigner distribution for parabolic waves was introduced earlier 8 and with it the correspond- 
ing radiative transfer equation has been derived with full mathematical rigor. 11 ' 12 In the 
case of anisotropic media fluctuating slowly in the longitudinal direction the 2f-RT equation 
developed here reduces to that of the paraxial waves in similar media which lends support to 
the validity of MSE. The other regime where the two frequency radiative transfer equation 
has been obtained with full mathematical rigor is geometrical optics. 13 

The main difference between the 2f-RT and the standard theory is that the former retains 
the wave nature of the process and is not just about energy transport. Hence the governing 
equation can not be derived simply based on the energy conservation law. 

2. Two-frequency Wigner distribution 

Let Uj,j = 1,2 be governed by the reduced wave equation 

AU J (r) + k*(v 1 + V 3 (r))U J (r) = f J (r), r G K 3 , j = l,2 (2) 

where Uj and Vj are respectively the mean and fluctuation of the refractive index associated 
with the wavenumber kj and are in general complex-valued. The source terms fj may result 
from the initial data or the external sources. Here and below the vacuum phase speed is set 
to be unity. To solve (T2J) one needs also some boundary condition which is assumed to be 
vanishing at the far field. 

We define the two-frequency Wigner distribution as 

In view of the definition, we see that both x and p are dimensionless. Here the choice of the 
scaling factors is crucial; namely, the spatial dependence of the wave field should be measured 
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w.r.t. the probing wavelength. The benefit is that this choice leads to a closed form equation 
for W . It is easy to see that the ensemble average (W) is just the (partial) Fourier transform 
of the mutual coherence function ([I]). The two- frequency Wigner distribution defined here 
has a different scaling factor from the one introduced for the parabolic waves. 8 

The purpose of introducing the two-frequency Wigner distribution is to develop a two- 
frequency theory in analogy to the well studied standard theory of radiative transfer. Al- 
though the definition Q requires the domain to be K 3 , the governing radiative transfer 
equation, once obtained, can be (inverse) Fourier transformed back to get the governing 
equation for the two-point function Lq(ri)[/J(r2) or ^12 as their boundary conditions are 
usually easier to describe (cf. eq. ( H2l) ). 

The Wigner distribution has the following easy-to-check properties: 

|2/„ „\J„J„ ( Vhh\ f {TT ,2/ X , f |rr 12/ 



|jy| 2 (x,p)dxdp = [^l^j J P\\ (x)dx J \U 2 \ 2 (x)dx 
I y(x,p)e^p = u^y-^-y.) (4) 
»'(x,p )e — <ix = ^ klk2 fu^ + ^f)Ui(^f- k -f), (5) 



where ^ stands for the Fourier transform, and hence contains all the information in the 
two-point two-frequency function. In particular, 

which, in the case of k\ = k 2 , is proportional to the energy flux density. 

We now derive the equation for the two-frequency Wigner distribution. After taking the 
derivative p • V and some calculation we have 

+ '-{ Vl -vl)W + F (6) 
where the function F depends linearly on Uj and ff 

+^ f e^U^ + ^-)/ 2 *(^ - i-)dy. (7) 



2(2tt) 3 j 

Substituting the spectral representation of Vj 



VG-(x) = / e^fydq) (8) 
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in the expression and using the definition of W we then obtain the exact equation 

p-VW- % -{v x -v* 2 )W-F (9) 
= l -j V 1 (dq)e^/ fc W(x,p - ±)-\J ^q)e- iqx / fc W(x,p - ^-). 

Here and below V^* is the complex-conjugate of the Fourier spectral measure V 2 - The full 
derivation of is given in Appendix A. 

Let us pause to compare the classical wave with the quantum wave function in the context 
of two-frequency formulation. The quantum wave functions \Pj at two different frequencies 
uji,uj 2 satisfy the stationary Schrodigner equaiton 

yA^+fe + ^x))^ = -Ujf&j + fj, J = 1,2, (10) 

where Uj + Vj are hypothetical, energy-dependent real-valued potentials. Here the source 
terms fj equal the initial data / of the time dependent problem. Usually in the quantum 
mechanical context, the potential function does not explicitly depend on the energy level 
(i.e. dispersionless). 

The natural definition of the two-frequency Wigner distribution for the quantum wave 
functions is 

^(x,p) = ^ J e-^ 1 (x+^)^(x-^)dy (11) 
which satisfies the Wigner-Moyal equation 

p- VW + i{u 2 - ujW + -vJW (12) 



h 



J Vi(dq)e*»' x ^(x, p - ^) - % - J V 2 *(dq)e-^W(x, p - ^) + F 



where F has a similar expression to ([?]). The main difference between the quantum and 
classical waves in the Wigner formulation is that the derivation of a closed-form equation 
does not require rescaling each energy component w.r.t. its de Broglie wavelength. The 
implication in radiative transfer will be further discussed (see the remark following eq. 027p ). 

3. Two-frequency radiative transfer scaling 

We assume that Vj(x),j = 1,2 are real-valued, centered, random stationary (i.e. statisti- 
cally homogeneous) ergodic field admitting the spectral representation (jSJ) with the spectral 
measures Vj(dp),j = 1,2 such that 



V 3 (dp)VUdq)) = 5(p-q)$ 3 -(p)dpdq 
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where $j are the (no nnegative- valued) power spectral densities of the random fields Vj, j = 
1,2. The above 5 function is a consequence of the statistical homogeneity of the random field 
Vj. As Vj, j — 1, 2 are real- valued, V*(dp) = Vj(— dp) and hence the power spectral densities 
$j(p) satisfy the symmetry property $j(p) = $.,•(— p), Vp. 

We will also need the cross-frequency correlation and we postulate the existence of the 
cross-frequency spectrum $ 12 such that 

(Vi(dp)t£(dq)) = 5(p - q)$ 12 (p)rfprfq. 

Here $12 needs not be real-valued. 

An important regime of multiple scattering of classical waves takes place when the scale 
of medium fluctuation is much smaller than the propagation distance but is comparable or 
much larger than the wavelength. 17 ' 19 Radiative transfer regime can be characterized by the 
scaling limit which replaces Vj + Vj in eq. (J2J) with 

^Lfuj + ^iVjC-)), ^>0, e<l (13) 

where e is the ratio of the scale of medium fluctuation to the 0(1) propagation distance and 
9 the ratio of the wavelength to the scale of medium fluctuation. Hence 9e is the ratio of the 
wavelength to the propagation distance and the prefactor {9e)~ 2 arises from rescaling the 
wavenumber k — > k/(e8). This is so called the weak coupling (or disorder) limit in kinetic 
theory which prohibits the Anderson localization from happening. 25 Note that the resulting 
medium fluctuation e~ 3 ^ 2 Vj(r / e) converges to a spatial white-noise in three dimensions. 

Physically speaking the radiative transfer scaling belongs to the diffusive wave regime 
under the condition of a large dimensionless conductance g = N£ t /L, where it is the transport 
mean free path, L is the sample size in the direction of propagation and N = 2nA/\ 2 is 
the number of transverse modes, limited by the illuminated area A and the wavelength of 
radiation A. 2 ' 23 The dimensionless conductance g can be expressed as g = k£t/Fr with the 
inverse Fresnel number Fr = XL/A. With the scaling (|13j) . kit ~ Fr _1 ~ 9~ 1 e~ 1 and hence 
g ~ 9~ 2 £~ 2 ^> 1 for any finite 9 as e — > 0. 

Anticipating small-scale fluctuation due to ffTSl we modify the definition of the two- 
frequency Wigner distribution in the following way 

Eq. ([9]) now becomes 

p-VW-F = ^-Ay x -vl)W + -=CW (14) 
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where the operator C is defined by 

£W(x,p) = ±J V^dqy^W^p-^-)-^ J y;(d q ) e -fw(x,p-|i). 

To capture the cross-frequency correlation in the radiative transfer regime we also need to 
restrict the frequency difference range 

limfci = lim/co = k, 2 — 1 = 3 (15) 

e-o e^o ' e6k v ' 

where k, 3 > are independent of e and 0. Assuming the differentiability of the mean 
refractive index's dependence on the wavenumber we write 

where v' is independent of e, 6. 

4. Multi-scale expansion (MSE) 

To derive the radiative transfer equation for the two-frequency Wigner distribution we em- 
ploy MSE 1 ' 20 which begins with introducing the fast variable 

x = x/e 

and treating x as independent from the slow variable x. Consequently the derivative p ■ V 
consists of two terms 

P- V = p- V x + e-V V s . (17) 
Then MSE posits the following asymptotic expansion: 

W(x,p) = iy(x,x,p) + v ^^ 1 (x,x,p)+^ 2 (x,x,p)+0(5 3/2 ), x = x^ 1 (18) 

whose proper sense will be explained below. 

Substituting the ansatz into eq. (fill) and using (ITTI) we determine each term of ([TBI by 
equating terms of the same order of magnitude starting with the highest order e" 1 . 

The e~ 1 -order equation has one term: 

p . Vxl? = 

which can be solved by setting W = W"(x, p). Namely, to the leading order W is independent 
of the fast variable. Since the fast variable is due to medium fluctuation, this suggests that 
W is deterministic. 
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The next is the e 1 / 2 -order equation: 

p . Vx^i = CW. (19) 

We seek a solution that is stationary in x, square-integrable in p and has finite second 
moment. The solvability condition (Fredholm alternative) is that the right hand side, CW, 
satisfies / E[&*CW]dp = for any x-stationary, square-integrable field \l/(x, p) satisfying 
p ■ Vx$ = 0. The solvability condition is, however, not easy to enforce. Alternatively we 
consider the regularized equation 

eW{ + p • VxWf = CW (20) 

which is always solvable for e > and admits the solution 

W?(x,x,p) = 4 / ^iW Jl W(x,p-^) (21) 

20 J e + iq ■ p/ki 2k\ 



-4 / e ~ lX lh w ^v-w)- 

29 J e — «q • p/fc2 2fc 2 



In the jargons of asymptotic analysis, 1 \feW[ is called the first corrector. In order to control 
the first corrector, we choose W such that CW has zero mean. This is a necessary condition 
as we seek a x-stationary solution and consequently (p • VxWi) = p • V x (Wi) = 0. Needless 
to say, this condition is weaker than the solvability condition stated above and is satisfied 
for any deterministic W since both V\ and Vi have zero mean. 

Indeed, under the assumption of deterministic W, the resulting equation will be much 
simplified so we impose this property on W from now on. The fact that in the limit W is 
deterministic can be proved rigorously in the paraxial regime. 12 

Finally the 0(1) equation is 



2h 



P- V x ^ 2 (x,x,p) = - p .VjV( x , p )-iv>W + F+^Jv 1 (dq)e i ^W?(x,%p 

~2l J ^*(dq)e-^W?(x,x,p - ^) (22) 

which can be solved with regularization as in ( 1201) and yields the second corrector eW^- Again 
we impose on the right hand side of (l2"2l) the weaker condition of zero mean. Using ( J2TT) in 
( |22|) . taking the ensemble average and passing to the limit e — > we obtain the governing 
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equation for W: 
p ■ V X W (x, p) + iv'W - (F) 



^ / ^(^(p - q))7r5(|p| 2 - |q| 2 Mx, p) + ^ f dq 







26 A 



iPl 



|q| 



W(x,p) 



Jl / ^ 2 (|(P " q))^(| P | 2 - |q| 2 Mx, p) - || 



4fl 2 



rfq$ 12 (q)e* i - 



2fci' 



+^ | c /q$ 12 (q)e- ^ 1 -^ 1 ) 7 r5(^--(p-|^))iy(x,p 

z 



4# 2 



<iq 



fc 2 VP 2fci^ fci 2fc 2 ^ 



$i 2 (q)e 



$ a (fr(p-q)) 

rfq lpl 2 -lql 2 (X ' P) 
9q 6q 

2k[ ~ 2k 2 ' 
Bq 9q 
2k[ ~ 2k 2 ' 

a-qcfcr 1 -^ 1 )^/ ^q ^q 



where we have used the fact that in the sense of generalized function 

1 % 
lim — = 7i5(£) — - 

with the second term giving rise to the Cauchy principal value integral denoted by -f . From 
(J7|) we have the expression for (F) 



(F) 



2(2tt) 3 
% 
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which depends only on the mean fields {U\) , (U2), both assumed known throughout the 
paper. 

Putting all the terms together with the regularization we arrive at the following MSE 



W(x t p) = W (x, p) + V£Wf(x, x, p) + £ ^(x, x, p) 



(23) 



which satisfies 

(p-V - -^=C^W + iv'W -F (24) 
= (zV - l)y/eWl + v^P • V x Wf - + (zV - l)eW| + ep • V X W|. 



Unfortunately the right hand side of ( l24|) does not vanish in the strong L 2 -topology but only 
in the weak topology as in 



lime / dx. 



dp^(x ) - I p)^(p) 



0, e 1} 



(25) 



(see Appendix B). It is not clear at this point how to justify the preceding argument and 
construction of asymptotic solution with full mathematical rigor. Fortunately, in the regime 
of geometrical optics, the rigorous asymptotic result can be obtained by a probabilistic 
method 13 and is the same as derived by MSE (see Section [6]). Another regime for which the 
asymptotic result can be fully justified is paraxial waves which we will turn to in the next 
section. 

Due to the assumption 015|) and the assumed continuous dependence of the medium fluc- 
tuation on the frequency we have lim*^! = lim$ 2 — lim$ 12 = $. As a consequence, all the 
Cauchy principal value integrals cancel out. With some changes of variables the governing 
equation for W takes the much simplified form: 

p • V X W + iv'W - (F) (26) 
irk 3 f . ,. ,k 
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J rfq$(^(p-q))5(|p| 2 - |q| 2 ) jMP-<i)P W ( x? q ) _ W ( x , p 



The 5-function in the scattering kernel is due to elastic scattering which preserve the 
wavenumber. When /3 = (then v\ = v-i and iv' ~ the imaginary part of u), eq. ( |26l) 
reduce to the standard form of radiative transfer equation for the phase space energy den- 
sity. 5 ' 16 ' 19 ' 22 For (3 > 0, the wave featue is retained in (j26|) . When (3 — > oo, the first term in 
the bracket on the right hand side of (|26|) drops out, due to rapid phase fluctuation, so the 
random scattering effect is pure damping: 

p-^W + iu'W-iF) = -^ydq$(^(p-q))5(|p| 2 -|q| 2 )W(x,p). 

As a comparison, for Schrodinger equation fTlTJT) in the frequency domain, we modify the 
Wigner distribution as 

^(x,p) = / e-^i(x+^)^(x-^)dy 

and in the limit e — > obtain the radiative transfer equation following the same procedure 

p-V^W + i(u 2 -uj 1 )W + ^u'W - (F) (27) 

n 

= ^ | rfq $(P T ^) ( 5(|p| 2 -| q | 2 )[ty(x,q)-iy(x,p 

The absence of the factor e lx '( p_q ^ in eq. ( |27l) . and therefore the cross-frequency interference, 
is the main characteristic of 2f-RT for quantum waves. 

5. Paraxial 2f-RT: anisotropic medium 

Forward-scattering approximation, also called paraxial approximation, is valid when back- 
scattering is negligible and, as we show now, this is the case for anisotropic media fluctuating 
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slowly in the (longitudinal) direction of propagation. Let z denote the longitudinal coordinate 
and x_i_ the transverse coordinates. Let p and denote the longitudinal and transverse 
components of p G IR 3 , respectively. Let q = (q, q_i_) G M 3 be likewise defined. 

Consider now a highly anisotropic spectral density for a medium fluctuating much more 
slowly in the longitudinal direction, i.e. replacing $((p — q)k/9) in ff26]) by 

^(^(p-g),^(p±-q±)), /7«1, 
which, in the limit rj — ► 0, tends to 

^(P -Q)J dw$ (w, ~(p ± - q ± )J . (28) 

Writing = ^(z, xj_, p, pj_) we can approximate eq. ( 1261) by 

= J dq ± J dw$(w,^(p ± - q±))5(\p±\ 2 - |q±| 2 ) 

e ^-( P X-qx)^(^ x±) ^ q± ) _^ )x±;p;p± )j. (29) 

Eq. ( 1291) is identical to the 2f-RT equation rigorously derived directly from the paraxial 
wave equation for similar anisotropic media. 11 ' 12 This is somewhat surprising in view of the 
different scaling factors in the definition of two-frequency Wigner distributions in the two 

cases. 

Note that in eq. (1291) the longitudinal momentum p plays the role of a parameter and does 
not change during propagation and scattering. An important implication of this observation 
is that eq. (|29|) can be solved as an evolution equation in the direction of increasing z with 
the one-sided boundary condition (e.g. at z = const.). In other words, the influence from 
the other boundary vanishes as the longitudinal direction is infinitely long. The initial value 
problem of ( 1291) is much easier to solve than the boundary value problem of ( 1261) . 



6. Two-frequency geometrical radiative transfer (2f-GRT) 

Let us consider the further limit 9 1 when the wavelength is much shorter than the corre- 
lation length of the medium fluctuation. To this end, the following form is more convenient 
to work with 

p ■ V X W + iv'W - (F) (30) 
= ^/^*(q)*(q-(p-^^ 
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which is obtained from eq. (1261) after a change of variables. We expand the right hand side 
of ( 1301) in 9 and pass to the limit 9 — > to obtain 

p ■ V X W + iv'W -(F) = — (V p - z/3x) • D ■ (V p - z'px) W (31) 
with the (momentum) diffusion coefficient 

D(p) =n J $(q)5(p ■ q)q ® qdq. (32) 

The symmetry $(p) = $(— p) plays an explicit role here in rendering the right hand side of 
eq. (130p a second-order operator in the limit 9 — > 0. Eq. ( 13TI) can be rigorously derived from 
geometrical optics by a probabilistic method. 13 

6. A. Spatial (frequency) spread and coherence bandwidth 

Through dimensional analysis, eq. (131]) yields qualitative information about important phys- 
ical parameters of the stochastic medium. To show this, let us assume for simplicity the 
isotropy of the medium, i.e. $(p) = $(|p|), so that D = C|p|~ 1 n(p) where 

C=l 1 5(^-AW|q|)|qMq (33) 



3J Mp| |q 

is a constant and II(p) the orthogonal projection onto the plane perpendicular to p. In 
view of (13"Tj) C (and D) has the dimension of inverse length while the variables x and p are 
dimensionless. 

Now consider the following change of variables 

x = a x kx, p = o-pp/k, (3 = fie/3 (34) 

where a x and a p are respectively the spreads in position and spatial frequency, and (3 C is the 
coherence bandwidth. Let us substitute (|34|) into eq. ( 13~TT) and aim for the standard form 

p • V^W + iv'W — (F) = (Vp-0xj • |p| _1 n(p) (Vp-j^x)f. (35) 

The 1-st term on the left side yields the first duality relation 

a x /a p ~ 1/k 2 . (36) 

The balance of terms in each pair of parentheses yields the second duality relation 

O'xO'p ~ -5- (37) 

whose left hand side is the space- spread-bandwidth product. Finally the removal of the con- 
stant C determines 

op ~ k 2/3 C 1/3 (38) 
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from which a x and f3 c can be determined by using (|36|) and (1371) : 

a x ~ AT^C 1 ^ & ~ A; 2 / 3 C- 2 / 3 . 

We do not know if, as it stands, eq. (1351) is analytically solvable but we can solve analytically 
for its boundary layer behavior. 

6.B. Boundary layer asymptotics: paraxial 2f-GRT 

Consider the half space z > occupied by the random medium and a collimated narrow- 
band beam propagating in the z direction and incident normal to the boundary (z = 0) of 
the medium. Near the point of incidence on the boundary the corresponding two-frequency 
Wigner distribution would be highly concentrated at the longitudinal momentum, say, p = 1. 
Hence we can assume that the projection n(p) in (135]) is effectively just the projection onto 
the transverse plane coordinated by xi and approximate eq. ( |3TT) by 

d z + p ± -V x± ]w + iv'W-{F) = ^(V p± -i(3x ± ) 2 W (39) 

where the constant C± is the paraxial approximation of (J32l) for \p\ — 1: 

C ± = ^ J $(0,q ± )|qj 2 dq ± . 

Here we have assumed the isotropy of $ in the transverse dimensions. Note that the longitu- 
dinal (momentum) diffusion vanishes and that the longitudinal momentum p plays the role 
of a parameter in eq. ( |39l) which then can be solved in the direction of increasing z as an 
evolution equation with initial data given at a fixed z. This is another instance of paraxial 
approximation. 

Let cr* be the spatial spread in the transverse coordinates xj_, i c the coherence length in the 
transverse dimensions and (3 C the coherence bandwidth. Let L be the scale of the boundary 
layer. We then seek the following change of variables 

x± = ^r, P± = p±k£ c , z = YT, P = 1T ( 40 ) 
a*/c Lk p c 



to remove all the physical parameters from (1391) and to aim for the form 

dgW + p± • V i± W + Lkiu'W - Lk (F) = (Vp ± - z/3x ± ) 2 W. (41) 
The same reasoning as above now leads to 

4<r* ~ L/k, a*/£ c ~ 1//3 C , 4 ~ k^L-^Cl 1 ' 2 



and hence 
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The layer thickness L may be determined by £ c ~ 1, i.e. L ~ k 2 C ± 1 . 
After the inverse Fourier transform eq. (jUJ) becomes 



w 9± • Vx ± r 



Lkiv'T - Lk (F) 



\y± + P*±\r 



(42) 



which is the governing equation for the two-frequency mutual coherence in the normalized 
variables. With data given on z = and vanishing far-field boundary condition in the 
transverse directions, Eq. (I42p can be solved analytically and its Green function is given by 

e -iLkv> (£40)1/2 r 1 



— exp 
(27r) 2 5sinh [{iA0)VH] 

coth [(iApy^z] 
(i4/3)V2 

tanh [(i4/3) 1/2 5] 



x exp 



ViA(5z 



y± - /3x ± -y' ± + f3x[ 



(43) 



cosh [(i4py/*z] 



x exp 



(z4/3)V2 



Formula (143|) is consistent with the asymptotic result in the literature which mainly con- 
cerns with the cross-frequency correlation of intensity. In the radiative transfer regime con- 
sidered here, the cross-spectral correlation of intensity is the square of the two-frequency 
mutual coherence and has the commonly accepted form 15 ' 24,26 



exp 



(44) 



which is just the large (3 asymptotic of the squared factor | sinh [(i4 ( 5) 1 ' 2 ^] |~ 2 in (143!) at 
z = 1 (see 13 for detailed comparison). Moreover (143]) provides detailed information about 
the simultaneous dependence of the mutual coherence on the frequency difference and spatial 
displacement for z e (0, l). 23 > 26 

Surprisingly, a closely related equation arises in the two-frequency formulation of the 
Markovian approximation of the paraxial waves. 8 The closed form solution is crucial for 
analyzing the performance of time reversal communication with broadband signals. 10 The 
solution procedure for (1431) is similar to that given elsewhere 10 and is omitted here. 



6.C. Paraxial 2/-GRT in anisotropic media 

We use here the setting and notation defined in Section [5] for anisotropic media. For simplicity 
we will set p — 1 and omit writing it out in W. In view of (T2SJ) we replace $(q) in (13"2l by 



5(q) J dw$(u7,q ± ) 
and obtain the transverse diffusion coefficient 

D_l(p_l) = 7T J dq± J dw§(w,q_L)5(p± ■ q±)q± <g> q± 
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whereas the longitudinal diffusion coefficient is zero. 

For simplicity we assume the isotropy in the transverse dimensions, Q(w, p_i_) = <&(w, |p_i_|), 
so that Dj_ = Cj_|p_i_| _1 Il_|_(p_|_) where 

C ± = £ / S(^j ■ t^tW, \q x \)\q ± \dwdq x 

1 j v ip±i |q±i J 

is a constant and IIx(px) is the orthogonal projection onto the transverse line perpendicular 
to p_i_. Hence eq. (1311) reduces to 



d z + p±- V x± W + ii/W - (F) (45) 

= ^ (V p± - i/?x_0 • IpxI^nxCpx) (V p± - i/3xx) 

Alternatively, eq. (1431) can also be derived from eq. (1291) by taking the geometrical optics 
limit as described in the beginning of Section 6. 

Consider the change of variables (1401) to remove all the physical parameters from (1431) and 
to aim for the form 

d s + px • V ix j + Lkiu'W - Lk (F) (46) 
= (V p± - ^xx) ■ |px|- 1 n ± ( P± ) (V px - z/3xx) IV 

where L should be interpreted as the distance of propagation. 
Following the same line of reasoning, we obtain that 

4a, ~ L/k, aj£ c ~ 1//3 C , 4 ~ C^L-^k^ 

and hence 

^~Cl /3 L 4 / 3 , /3 C ~ C- 2/3 L- 5 / s k-\ 
Unlike (139|) it is unclear if a closed- form solution to eq. (j4"5l) exists or not. 

7. Discussion and conclusion 

The standard (one-frequency) RT can be formally derived from the wave equation in at 
least two ways: the diagrammatic expansion method, as the ladder approximation of the 
Bethe-Salpeter equation, 19 ' 26 and the multi-scale expansion method advocated here. 1 The 
latter is considerably simpler than the former in terms of the amount of calculation involved. 
Both approaches have been developed with full mathematical rigor in some special cases 
(see 7,9 and the references therein). There are two regimes for which the 2f-RT equation has 
been derived with full mathematical rigor: first, for the paraxial wave equation by using the 
so called martingale method in probability theory; 11,12 second, for the spherical waves in 



11 



geometrical optics by the path-integration method. 13 These rigorous results coincide with 
those derived here for the respective regimes and hence support the validity of MSE. 

Within the framework of 2f-RT, a paraxial form arises naturally in anisotropic media 
which fluctuate slowly in the longitudinal direction. Another form of paraxial 2f-RT takes 
place in the boundary layer asymptotics of isotropic media. The latter equation turns out 
to be exactly solvable and the boundary layer behavior is given in a closed form, revealing 
highly non-trivial structure of the two-frequency mutual coherence. In any case, dimensional 
analysis with the 2f-GRT equations yields qualitative scaling behavior of the spatial spread, 
the spatial frequency spread and the coherent bandwidth in various regimes. 

From the point of view of computation, especially Monte Carlo simulation, it appears to 
be natural to introduce the new quantity 

2U(x, p) = e~ mp W(x, p) 

and rewrite eq. ( 1261) in the following form 

p ■ V X 2U + if3\p\ 2 W + iu'W - e" i/3x - p (F) 

2U(x,q) -2U(x,p) . 

The solution 2U can then be expressed as a path integration over the Markov process gener- 
ated by the operator A defined by 

2U(x,q) -2U(x,p) 

when V is real-valued and $ is nonnegative. We will pursue this observation in a separate 
publication. 13 
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A. Derivation of eq. ([6]) 

Applying the operator p ■ V to the definition ([3]) we obtain 
Integrating by parts with the first V y in the above integrals, we have 



P-V X W = 




where the other resulting terms are canceled with each other. From eq. (J2J), 




Using (H8~j) in (HT) we arrive at eq. (JBj). 
B. Weak convergence of corrector 

First we show that the corrector does not vanish in the in the mean-square norm in any 
dimension, i.e. lim e _>o£ / (|Wf | 2 ) rfxdp > in general. For simplicity, consider only the term 
involving, say, V± in the expression (T2TT) . A straightforward calculation shows 

l[m nJw I dpdxdq^iq) £ |iy( x ,p-^)| 2 

= ^ J rfpdxrfq$(q)(5(p-q)|iy(x,p- ^)f 
which is positive in general. 
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Next we show that the corrector vanishes in the weak topology 



lime / dx. 

£->0 



dp^ 1 e (x, 7 ,p)V(p) 



0, G L 2 



(49) 



It suffices to prove (149j) for any smooth, compactly supported function ^. For the term 
involving Vi only, we have 



$i(qMp)V*(p') 



lim — — / dpdp' (ixrfq . 

e^oAO 2 ] yy ^(£ + z'p-q/A;i)(e-ip / -q/A;i 

,2 



*o 4# 2 



vr / dp5(p-q)^(p)^(x,p-^) 



dp x,p - — 

p • q 2ki 



vr / dp'5( p '.q)^(p')^(x,p'-^) 



where q = q/|q| for sufficiently smooth W, $ and rapidly decaying $. The essential point now 
is that |q| -2 is an integrable singularity in three dimensions and hence the above expression 
vanishes in the limit. 
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